/*
 *  smoother.cpp
 *
 *
 *  Created by Steve Fan on 13/08/08.
 *
 */
#include <math.h>
#include <stdio.h>
#include <S.h>

#define THREE_PI (3*M_1_PI)
#define R64

#ifdef R64
typedef int intType;
#else
typedef long intType;
#endif

inline double ifx(double y, double lb, double ub);
inline double igy(double x, double lb, double ub);
inline double ifgxy(double lb, double ub);
double condA(double x, double y);
double condB(double x, double y);
double condC(double x, double y);
double condD(double x, double y);
double condE(double x, double y);
double condF(double x, double y);
double condG(double x, double y);
double condH(double x, double y);
double condI(double x, double y);
double pbivarbiwgt(double x, double y);


extern "C" void smoother_exact(double *cx, double *cy, intType *ncell, double *quadPts, double *quadWgt,
	intType *nquadPts, double *bandwidth, double *cellwidth, double *offsets, double *Lambda,
	double *xobs, double *yobs, intType *nobs) {

	intType rowCell, colCell, R, C, i; //ievalpt;
	double x, y, xdist, ydist, dist, xleft, ybottom, xright, ytop, denominator, tmp_integral, numerator;
	double maxdist = *bandwidth + 2*M_SQRT2*(*cellwidth), halfcw = (*cellwidth)/2;

	for (rowCell=0; rowCell<*ncell; rowCell++) {
		for (C=0; C<*nquadPts; C++) {
			for (R=0; R<*nquadPts; R++) {
				numerator = 0;
				denominator = 0;
				x = cx[rowCell] + halfcw*quadPts[R];
				y = cy[rowCell] + halfcw*quadPts[C];

        			for (i=0; i<*nobs; i++) {
          				xdist = (xobs[i] - x)/(*bandwidth);
          				if (fabs(xdist) > 1) ;
          				else {
						ydist = (yobs[i] - y)/(*bandwidth);
          					if (fabs(ydist) > 1) ;
						else {
							dist = sqrt(pow(xdist, 2) + pow(ydist, 2));
							if (dist > 1) ;
							else {
                						numerator+= (pow( (1- pow(dist, 2)), 2)*THREE_PI);
							}
            					}
					}
				}
        
				for (colCell=0; colCell<*ncell; colCell++) {
					xdist = cx[colCell] - x;
					if (fabs(xdist) > maxdist) ;
					else {
						ydist = cy[colCell] - y;
						if (fabs(ydist) > maxdist) ;
						else {
							dist = sqrt(pow(xdist, 2) + pow(ydist, 2));
							if (dist > maxdist) ;
							else {
								xleft    = (xdist - halfcw)/(*bandwidth);
								ybottom  = (ydist - halfcw)/(*bandwidth);
								xright	 = (xdist + halfcw)/(*bandwidth);
								ytop     = (ydist + halfcw)/(*bandwidth);
								tmp_integral = (pbivarbiwgt(xright, ytop) + pbivarbiwgt(xleft, ybottom)) -
									(pbivarbiwgt(xleft, ytop) + pbivarbiwgt(xright, ybottom));
								if ((tmp_integral > 0) && (offsets[colCell] > 0)) denominator += (offsets[colCell]*tmp_integral);
							}
						}
					}
				}

				Lambda[rowCell] += (quadWgt[R]*quadWgt[C]/4)*(numerator/denominator);
			}
		}
	}
}

inline double ifx(double y, double lb, double ub) {
	// \int_{lb}^{ub} f(x, y) \, \mathrm{d}x
	return(-1*((pow(y, 5)/5 - 2*pow(y, 3)/3 + y/2)*(ub-lb) + (pow(ub,3)-pow(lb,3))*pow(y, 3)/9));
}

inline double igy(double x, double lb, double ub) {
	// \int_{lb}^{ub} g(x, y) \, \mathrm{d}y
	return((pow(x,5)/5 - 2*pow(x, 3)/3 + x/2)*(ub-lb) + (pow(ub,3)-pow(lb,3))*pow(x, 3)/9);
}

inline double ifgxy(double lb, double ub) {
	// \int_{\dai C} f(x, y) \, \mathrm{d}x + g(x, y) \, \mathrm{d}y in term of polar
	// coordinate with r^2 = 1
	// lb is the lower angle; ub is the upper angle
	return((5*(ub-lb) -( sin(4*ub)-sin(4*lb) ))/30);
}

double condA(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), y) +
            ifx(y, x, -1*sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condB(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), sqrt(1-pow(x, 2)))+
            ifgxy(acos(x), asin(y)) +
            ifx(y, sqrt(1-pow(y, 2)), -1*sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condC(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), sqrt(1-pow(x, 2))) +
            ifgxy(acos(x), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condD(double x, double y) {
  double result;
  result = ifx(y, sqrt(1-pow(y, 2)), -1*sqrt(1-pow(y, 2))) +
              ifgxy(M_PI-asin(y), 2*M_PI+asin(y));
  return(THREE_PI*result);
}

double condE(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), y) +
            ifx(y, x, -1*sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condF(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), y) +
            ifx(y, x, -1*sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condG(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), sqrt(1-pow(x, 2))) +
            ifgxy(acos(x), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condH(double x, double y) {
  double result;
  result = igy(x, -1*sqrt(1-pow(x, 2)), y) +
            ifx(y, x, -1*sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI-acos(x));
  return(THREE_PI*result);
}

double condI(double x, double y) {
  double result;
  result = ifx(y, sqrt(1-pow(y, 2)), -sqrt(1-pow(y, 2))) +
            ifgxy(M_PI-asin(y), 2*M_PI+asin(y));
  return(THREE_PI*result);
}

double pbivarbiwgt(double x, double y) {
	double xsq = pow(x, 2), sqrtonemxsq = sqrt(1-xsq), ans;

	if (x<= -1)                           		{ans = 0;}
	else if ((-1<x) && (x<=0)) {
    	if (y <= -sqrtonemxsq)              	{ans = 0;}
    	else if ((-sqrtonemxsq<y) && (y<=0)) 	{ans = condE(x, y);}
    	else if ((0<y) && (y<=sqrtonemxsq))  	{ans = condF(x, y);}
    	else                                	{ans = condG(x, y);} // (y>sqrt(1-x^2))
	}
	else if ((0<x) && (x<=1)) {
    	if (y <= -1)                        	{ans = 0;}
    	else if ((-1<y) && (y<=-sqrtonemxsq))  	{ans = condI(x, y);}
    	else if ((-sqrtonemxsq<y) && (y<=0)) 	{ans = condH(x, y);}
    	else if ((0<y) && (y<=sqrtonemxsq))  	{ans = condA(x, y);}
    	else if ((sqrtonemxsq<y) && (y<=1))  	{ans = condB(x, y);}
    	else                                	{ans = condC(x, y);}// (y>1)
	}
	else { // x>1
		if (y<=-1)                           	{ans = 0;}
		else if ((-1<y) && (y<=0))            	{ans = condI(x, y);}
		else if ((0<y) && (y<=1))             	{ans = condD(x, y);}
		else                                 	{ans = 1;} // y>1
	}
	return(ans);
}

